一种基于压缩感知的低空风切变风速估计方法。其包括以下步骤:1)利用多普勒矢量构建冗余字典,以实现回波信号的稀疏表示;2)构造测量矩阵;3)计算信号压缩后的测量值;4)利用正则化正交匹配追踪算法重构稀疏信号;5)利用复幅度估计值中的非零值获得雷达回波信号的多普勒频率估计;6)计算杂波抑制后复幅度估计值中最大值的位置,即为风切变信号的风速估计结果;7)重复步骤3)-6),依据距离单元数判断风速估计是否完成,依次对范围内的所有距离单元的回波数据进行风速估计。本发明针对脉冲数较少且信噪比较低时风速估计精度变差的问题,提出了一种基于压缩感知的低空风切变风速估计方法。该发明方法能够在获得精确风速估计的同时使得频谱分辨率大大提高,即能够很好地区分在频域间隔很近的风切变与地杂波信号。
1.一种基于压缩感知的低空风切变风速估计方法,其特征在于:针对某一距离单元的回波数据,所述的风速估计方法包括按顺序进行的下列步骤:1)根据机载气象雷达回波信号的稀疏性,利用多普勒矢量构建冗余字典,以实现回波信号的稀疏表示;2)构造测量矩阵;3)计算信号压缩后的测量值;4)利用正则化正交匹配追踪算法重构稀疏信号;5)利用复幅度估计值中的非零值获得雷达回波信号的多普勒频率估计;6)计算杂波抑制后复幅度估计值中最大值的位置,即为风切变信号的风速估计结果;7)重复步骤3)—6),依据距离单元数判断风速估计是否完成,依次对范围内的所有距离单元的回波数据进行风速估计;所述的步骤4)中的正则化正交匹配追踪算法的实现步骤主要包括:利用相关性筛选原子,计算相关系数中K个最大值所对应的索引值,将其存入候选集中;将候选集中索引值对应原子的相关系数拆分成若干组进行正则化处理,选择能量最大的一组,即实现利用正则化原则筛选信号支撑点;更新支撑集和索引值集合;采用最小二乘法进行信号逼近,并更新余量;当索引值集合中的原子数大于等于两倍的稀疏度K时,停止迭代,即实现信号重构。
2.根据权利要求1所述的基于压缩感知的低空风切变风速估计方法,其特征在于:所述的步骤1)中对回波信号进行稀疏表示时,其稀疏度为K,即回波信号在稀疏域的表示系数只有K个非零值,并且K远小于信号长度。
3.根据权利要求1所述的基于压缩感知的低空风切变风速估计方法,其特征在于:所述的步骤2)中的测量矩阵为高斯随机矩阵。
4.根据权利要求1所述的基于压缩感知的低空风切变风速估计方法,其特征在于:所述的步骤5)中多普勒频率估计包括风切变信号和地杂波信号的多普勒频率估计。
5.根据权利要求1所述的基于压缩感知的低空风切变风速估计方法,其特征在于:所述的步骤6)中杂波抑制的方法是:利用回波信号的多普勒频率分布特性,将位于零多普勒频率附近的估计值置零。
6.根据权利要求1所述的基于压缩感知的低空风切变风速估计方法,其特征在于:所述的步骤7)中依次对范围内的所有距离单元的回波数据进行风速估计的方法是:考虑波束照射区域所对应的最大距离单元以内的雷达回波数据,从最近的距离单元开始计算。技术领域
本发明属于机载气象雷达低空风切变检测技术领域,特别是涉及一种基于压缩感知的低空风切变风速估计方法。
背景技术
作为对航空运输安全威胁最大的天气现象之一,低空风切变主要是指在航空器起降阶段,在相同或不同高度(在600m左右)的较短距离内,风向或风速发生快速变化的天气现象。低空风切变发生突然并且具有持续时间短、影响强度大、危险性高等特点。由于在飞机起降阶段可操纵的高度空间受限,若风切变强度较大就会极易引发严重的坠机事故。因此,低空风切变探测技术的研究已成为现代航空运输领域的一个重要课题。2004年,由Donoho、Candes和Tao等人提出的压缩感知(Compressivexa0Sensing,CS)理论是一种充分利用信号的稀疏性或可压缩性来实现信号处理的基础理论。该理论的应用提供了降低信号采样率的可能性,实现了信号采样理论的深刻变革,并具有广阔的应用前景,目前其研究成果已应用在光学成像、模拟信息转换、生物传感、图像融合、SAR目标识别、波形信号仿真分析、探地雷达成像、雷达目标探测和目标波达方向(Directionxa0ofxa0Arrival,DOA)估计等诸多领域。机载气象雷达下视探测低空风切变时,接收到的回波数据无可避免地要受到地杂波信号的影响。一般来说,由于相对运动速度的不同,回波信号的频谱中包括风切变信号和地杂波信号产生的两个频率谱峰,简称为双峰谱。依据风速估计时是否采用杂波抑制滤波器,常用的风速估计方法分为两大类。当功率谱为双峰谱时,即考虑在不抑制杂波情况下的信号频谱,常采用的方法包括基于参数化模型的谱估计方法和利用模式分析的扩展Prony方法;当功率谱为单峰谱时,即考虑没有杂波影响或抑制杂波后的信号频谱,应用最广的估计方法是在时域分析基础上利用相关信息的脉冲对(Pulsexa0Pairxa0Processing,PPP)法和基于频域分析的快速傅里叶变换(Fastxa0Fourierxa0Transform,FFT)法。总体来说,由于算法简单,运算量及硬件实现的设备量较小等优点,因此PPP法在气象雷达信号处理中的应用较广,但当信噪比较低时,其估计性能变差。由于FFT法在频域进行分析,因此能够较好地分辨非气象回波,在地杂波抑制方面要优于脉冲对法。但是,由于FFT法的谱分辨力受数据长度的限制,特别是在数据长度短和低信噪比时,FFT法和PPP法的风速估计性能均会变差。机载气象雷达回波信号中三个重要的气象目标参数是:信号功率、平均速度和谱宽。其中,平均速度反映了气象目标的运动特征。在实际的机载气象雷达风切变检测系统中,水平风速可由雷达回波的多普勒频移来获取。针对脉冲数较少且信噪比较低时,现有的FFT法和PPP法等风速估计方法存在谱分辨力受限,估计性能下降等技术难题。因此,研究在脉冲数较少且信噪比较低时的风速估计问题具有十分重要的现实意义。
发明内容
为了解决上述问题,本发明的目的在于提供一种在脉冲数较少且信噪比较低时,仍然具有较高风速估计精度的基于压缩感知的低空风切变风速估计方法。为了达到上述目的,本发明提供的基于压缩感知的低空风切变风速估计方法包括按顺序进行的下列步骤:1)根据机载气象雷达回波信号的稀疏性,利用多普勒矢量构建冗余字典,以实现回波信号的稀疏表示;2)构造测量矩阵;3)计算信号压缩后的测量值;4)利用正则化正交匹配追踪算法,重构稀疏信号;5)利用复幅度估计值中的非零值,获得雷达回波信号的多普勒频率估计;6)计算杂波抑制后复幅度估计值中最大值的位置,即为风切变信号的风速估计结果;7)重复步骤3)-6),依据距离单元数判断风速估计是否完成,依次对范围内的所有距离单元的回波数据进行风速估计。所述的步骤1)中对回波信号进行稀疏表示时,其稀疏度为K,即回波信号在稀疏域的表示系数只有K个非零值,并且K远小于信号长度。所述的步骤2)中的测量矩阵为高斯随机矩阵。所述的步骤4)中的正则化正交匹配追踪算法的实现步骤主要包括:利用相关性筛选原子,计算相关系数中K个最大值所对应的索引值,将其存入候选集中;将候选集中索引值对应原子的相关系数拆分成若干组进行正则化处理,选择能量最大的一组,即实现利用正则化原则筛选信号支撑点;更新支撑集和索引值集合;采用最小二乘法进行信号逼近,并更新余量;当索引值集合中的原子数大于等于两倍的稀疏度K时,停止迭代,即实现信号重构。所述的步骤5)中多普勒频率估计包括风切变信号和地杂波信号的多普勒频率估计。所述的步骤6)中地杂波抑制的方法是:利用回波信号的多普勒频率分布特性,将位于零多普勒频率附近的估计值置零。所述的步骤7)中依次对范围内的所有距离单元的回波数据进行风速估计的方法是:考虑波束照射区域所对应的最大距离单元以内的雷达回波数据,从最近的距离单元开始计算。本发明提供的基于压缩感知的低空风切变风速估计方法具有如下优点和积极效果:与现有技术相比,本发明方法根据雷达回波中风切变信号相对于背景的稀疏性,利用多普勒矢量构建一个冗余字典以实现信号的稀疏表示,采用观测矩阵对信号进行压缩处理,并通过信号重构算法恢复该稀疏信号,实现风速的精确估计。当脉冲数较少且信噪比较低时,本发明方法能够在获得精确风速估计的同时使得频谱分辨率大大提高,即能够很好地区分在频域间隔很近的风切变与地杂波信号。
附图说明
图1为压缩感知理论的基本流程图。图2为本发明提供的基于压缩感知的低空风切变风速估计方法流程图。图3为FFT方法与本发明方法估计信号的频谱对比图。图4为利用本发明方法抑制杂波后的回波信号频谱分布图。图5(a)为脉冲数为50时全距离单元的风速估计结果。图5(b)为脉冲数为50时风切变信号风速估计局部图。图6(a)为脉冲数为32时全距离单元的风速估计结果。图6(b)为脉冲数为32时风切变信号风速估计局部图。
具体实施方式
下面结合附图和具体实施例对本发明提供的基于压缩感知的低空风切变风速估计方法进行详细说明。如图2所示,本发明提供的基于压缩感知的低空风切变风速估计方法包括按顺序进行的下列步骤:1)根据机载气象雷达回波信号的稀疏性,利用多普勒矢量构建冗余字典,以实现回波信号的稀疏表示;假设检测低空风切变所用的相干处理脉冲数为N_r(即,信号长度为N_r),脉冲重复频率为F_r,将多普勒频率范围(0-F_r)等分成N个离散的频点f₁,f₂,...,f_N。其中,j为虚数单位,频点f_i的多普勒矢量可以表示为: a ( f i ) = exp { j 2 π F r [ 0,1 , . . . , N r - 1 ] T f i } - - - ( 1 ) ]]>由各频点处的多普勒矢量构成的矩阵可写成:ntent="drawing" img-format="tif" inline="no" orientation="portrait" wi="569"/>ntent="drawing" img-format="tif" inline="no" orientation="portrait" wi="700"/>低空风切变由大量的微粒组成,属于分布式的气象目标。考虑强地杂波环境下的目标回波信号,将每个点散射体近似等效为理想的球体,计算每个散射单体的雷达回波,通过相干叠加生成每个距离单元的回波信号。因此,整个波束扫描范围内所有散射点回波的叠加即为雷达回波,可表示为:ntent="drawing" img-format="tif" inline="no" orientation="portrait" wi="700"/>其中,A_i(t)表示第i个散射点的回波幅度;ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="93"/>表示第i个散射点的回波相位;c(t)表示地杂波信号,n(t)表示干扰噪声信号,j表示虚数单位。将回波数据中的风切变信号和地杂波信号等效为信号S(t),即S(t)=s(t)+c(t)。雷达回波信号的矩阵形式为:x=S+nxa0xa0xa0xa0xa0(4)由各频点处的多普勒矢量构成的矩阵A即为对信号进行压缩处理的冗余字典(也可称为过完备库),其中矩阵A的每一列都是过完备库中的原子。利用过完备库中的原子代替传统的正交基,并找到具有最佳线性组合的m项原子来表示目标信号,即可实现在过完备库下的信号稀疏表示。因此,接收到的目标信号S可表示为A中原子的线性组合,即:x=Aα+nxa0xa0xa0xa0xa0(5)其中,α表示信号在频率域的复幅度,且α中只有K个非零分量,即雷达回波信号的稀疏度为K,且K≤N_r,所以α是稀疏的,实现了接收数据的稀疏表示。2)构造一个M×N_r维的高斯随机矩阵,并对每一列进行归一化得到测量矩阵Φ;3)利用测量矩阵Φ计算信号压缩后的M×1维的测量值y:y=Φx=ΦAα+Φnxa0xa0xa0xa0xa0(6)=Ψα+w式中,Ψ=ΦA为等效的测量矩阵,w=Φn为等效的噪声分量。4)利用正则化正交匹配追踪算法重构稀疏信号;正则化正交匹配追踪算法的基本思想是:挑选多个原子作为候选集,并按照正则化原则从中挑选出部分原子,即实现利用正则化原则筛选信号支撑点。在此基础上,将挑选的原子索引值并入最终的支撑集以实现原子的快速、有效选择。其具体实现方法包括:a)初始化余量r₀=y,估计信号的稀疏度为K,索引值集合ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="157"/>候选集ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="143"/>迭代次数t=1;b)利用相关性筛选原子,并计算余量r与等效的测量矩阵Ψ中各个原子之间内积的绝对值,利用式(7)来计算相关系数u,寻找u中K个最大值所对应的索引值,将其存入候选集J中;u={u_j|u_j=|<r,Ψ_j>|,j=1,2,...N}xa0xa0xa0xa0xa0(7)c)将候选集J中索引值对应原子的相关系数拆分成若干组进行正则化处理,选择能量最大的一组,将其相关系数对应的原子索引值存入候选集的索引值集合J₀中,即实现利用正则化原则筛选信号支撑点,且该集合中原子的相关系数必须满足式(8);|u(i)|≤2|u(j)|,i,j∈Jxa0xa0xa0xa0xa0(8)d)更新支撑集Φ_Λ,更新索引值集合Λ=Λ∩J₀;e)采用最小二乘法进行信号逼近,并更新余量; α ^ = arg min | | y - Φ Λ α | | 2 ]]>xa0xa0xa0xa0xa0(9) r new = y - Φ Λ α ^ ]]>f)若索引值集合中的原子数|Λ|≥2K,则停止迭代;否则令r=r_new,t=t+1,转到步骤b)继续执行直到停止迭代。5)利用正则化正交匹配追踪算法实现信号重构后,复幅度估计值ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="30"/>中非零值的位置即为雷达回波信号(包括风切变信号和地杂波信号)的多普勒频率估计。6)信号重构后得到的复幅度估计值ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="54"/>只存在零值和非零值两种状态。其中,非零值的位置即对应风切变信号和地杂波信号的速度估计。由于风切变信号和地杂波信号的相对速度不同,接收的雷达回波信号中强杂波主要集中在零多普勒频率附近,风切变信号的多普勒频率沿距离呈现出反“S”特性,故将位于零多普勒频率附近的估计值置零,即可实现地杂波抑制。杂波抑制后,复幅度估计值ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="30"/>中最大的非零值的位置即为风切变信号的风速估计结果。7)重复步骤3)-6),依据距离单元数判断风速估计是否完成,依次对范围内的所有距离单元的回波数据进行风速估计,方法是考虑波束照射区域所对应的最大距离单元以内的雷达回波数据,从最近的距离单元开始计算。假设k表示当前的距离单元数,且1≤k≤k_r,K_r表示机载气象雷达回波的最大距离单元数。图1为压缩感知理论的基本流程图。如图1所示,压缩感知的基本思想是:如果某种变换之后信号是可稀疏分解或可压缩的,则信号可投影到一个与变换基不相关的测量矩阵上,并获得远少于信号长度的测量值,最终通过求解优化问题精确地重构信号。图2为本发明提供的基于压缩感知的低空风切变风速估计方法流程图。其中关键步骤在于:依据距离单元数判断风速估计是否完成,考虑波束照射区域所对应的最大距离单元以内的雷达回波数据。以某一距离单元的回波数据为例,本发明方法是首先利用多普勒矢量构建一个冗余字典,并用过完备原子库中的原子实现信号的稀疏表示;其次,选择高斯随机矩阵作为观测矩阵对信号进行压缩;再次,通过信号重构算法估计回波信号的多普勒频率;最后,在分析风切变信号频谱与地杂波频谱特点的基础上,通过将零多普勒频率附近的估计值置零的方法实现地杂波抑制,并获得精确的低空风切变风速估计。参照实际装机使用的WXR-2100型机载气象雷达,基本的仿真参数设定如下:飞机高度600m,飞机速度75m/s,雷达发射功率100W,波束宽度3.5°,波束下视角2.75°,脉冲重复频率7000Hz,发射脉宽1μs,天线增益30dB,相干脉冲数128,杂噪比30dB,信噪比5dB。图3为FFT方法与本发明方法估计信号的频谱对比图。设脉冲数减少时的信号为顺序选取原始信号的前64个脉冲回波,以第49号距离单元(约7.35km)的回波数据为例,由图3可见,蓝色虚线和红色实心点曲线分别表示FFT方法与本发明方法估计的信号频谱。当相干处理的脉冲数较少时,频域能量积累不足且频谱展宽严重,在风切变信号频谱与地杂波信号频谱相近的区域,FFT方法的估计结果已无法区分风切变信号,但本发明方法仍能够清楚地表征原始信号的多普勒中心频率的频谱特征,即可获得风切变信号与地杂波信号的双峰谱。将估计结果与原始的风切变信号频谱和原始的地杂波信号频谱对比,可以发现频谱中的双峰位置与原始信号一致,更加证实了本发明提供的基于压缩感知的多普勒频率估计方法的有效性。图4为利用本发明方法抑制杂波后的回波信号频谱分布图。在回波信号频谱特征分析的基础上,利用强杂波主要集中在零多普勒频率附近的频谱分布特性,将位于零多普勒频率附近的多普勒估计值置零,从而实现抑制地杂波的目的。本发明方法避免了复杂的滤波器设计环节,能够简单有效地实现杂波抑制。图5为脉冲数为50时的速度估计对比图。图5(a)为脉冲数为50时全距离单元的风速估计结果。由图中可见,大约在0-2km和11-15km处不存在风切变信号,即只包含噪声和杂波信号,可以认为杂波抑制后的估计值为噪声速度。图5(b)为脉冲数为50时风切变信号风速估计局部图。该图可有助于重点观察存在风切变信号区域(约3-11km处)的风速估计结果。图6为脉冲数为32时的速度估计对比图。图6(a)为脉冲数为32时全距离单元的风速估计结果。图6(b)为脉冲数为32时风切变信号风速估计局部图。首先,分别利用固定凹口滤波器和自适应频域置零等方法抑制地杂波;其次,采用PPP法对杂波抑制后的回波数据进行风速估计;最后,将该估计结果与利用本发明方法得到的风速估计结果进行对比分析,由仿真结果可见:(1)当脉冲数减少到N_r=50时,利用固定凹口滤波器和自适应频域置零方法抑制杂波后的速度估计精度均有明显下降。相比之下,本发明方法的速度估计结果与原始的速度分布最为接近,估计误差较小。(2)当脉冲数减少到N_r=32时,由于时域数据减少,频域能量积累不足且频谱展宽严重,自适应频域置零方法抑制杂波后的速度估计已完全失真。固定凹口滤波器方法在风切变信号与地杂波信号频谱邻近的区域,速度估计产生了较大的波动。此时,由于该类方法的谱宽估计值较大,造成滤波器凹口增大,使得在抑制地杂波时风切变信号的能量损失严重,从而产生较大的估计误差。本发明提供的基于压缩感知的低空风切变风速估计方法对相干处理的脉冲数不敏感,当脉冲数减少到N_r=32时,本发明方法仍能获得高精度的风速估计,明显优于上述两种方法。